/*
    <one line to give the program's name and a brief idea of what it does.>
    Copyright (C) <year>  <name of author>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#include "qcCfour.h"

namespace atomism
{
 

const string qcCfour::BASISSET_NAMES[] =
{	"6-31G",	
	"6-311G",
	"ANO0",
	"ANO1",
	"ANO2",
	"AUG-pVDZ",
        "AUG-pVTZ",
        "AUG-pVQZ",
	"AUG-pCVDZ",
        "AUG-pCVTZ",
        "AUG-pCVQZ"
};

const string qcCfour::CITYPE_NAMES[] =
{	"HF",	
	"MP2",
	"MP4",
	"CCSD",
	"CCSD(T)",
	"CCSDT"
 };


bool qcCfour::isCfourRegisteredInPython=0;
msRegistrar qcCfour::Registrar("qcCfour", qcCfour::createInstance);
 
void qcCfour::registryInPython()
{
 #if USE_PYTHON
 
 msTreeMapper::registryInPython();
 
 if( ! qcCfour::isCfourRegisteredInPython ) {
   
     using namespace boost::python;
     class_<qcCfour, bases<msTreeMapper>, boost::shared_ptr<qcCfour> >(
           "qcCfour", no_init)
          .def( "New" ,	&qcCfour::New )
          .staticmethod("New")
	  .def("singlePoint" ,  (double (qcCfour::*)(const msZMat&) ) &qcCfour::singlePoint , 
	       "realize a single point energy calculation and return the energy. arg2: z-matrix")
	  .def("setLevel" ,  &qcCfour::setLevel , 
	       "set the level of theory. arg2: basis set (enum BasisSet), arg3: CI level (enum CiType)")
	  .def("setGenBasSymLink",&qcCfour::setGenBasSymLink,
               "set the symbolic link path for the file GENBAS");
	   enum_<qcBasisSet>("BasisSet")
              .value( "G6_31", G6_31)
              .value( "G6_311", G6_311)
	      .value(BASISSET_NAMES[ANO0].c_str(), ANO0)
	      .value(BASISSET_NAMES[ANO1].c_str(), ANO1)
	      .value(BASISSET_NAMES[ANO2].c_str(), ANO2)
	      .value("AUG_PVDZ", AUG_PVDZ)
	      .value("AUG_PVTZ", AUG_PVTZ)
	      .value("AUG_PVQZ", AUG_PVQZ)
	      .value("AUG_PCVDZ", AUG_PCVDZ)
	      .value("AUG_PCVTZ", AUG_PCVTZ)
	      .value("AUG_PCVQZ", AUG_PCVQZ)
              .export_values();
	   enum_<qcCiType>("CiType")
              .value(CITYPE_NAMES[HF].c_str(), HF)
              .value(CITYPE_NAMES[MP2].c_str(), MP2)
	      .value(CITYPE_NAMES[MP4].c_str(), MP4)
	      .value(CITYPE_NAMES[CCSD].c_str(),CCSD)
	      .value("CCSDpT", CCSDpT)
	      .value(CITYPE_NAMES[CCSDT].c_str(), CCSDT)
              .export_values();
    qcCfour::isCfourRegisteredInPython=1;
    }
#endif
 }


//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

boost::shared_ptr<msTreeMapper> qcCfour::setLevel(qcBasisSet bs, qcCiType ci){

 BasisSet=BASISSET_NAMES[bs];CiType=CITYPE_NAMES[ci];
 return mySharedPtr();
 }

//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

boost::shared_ptr<msTreeMapper> qcCfour::setWorkingDirectory(std::string dir){

 WorkingDirectory=dir;
 return mySharedPtr();
 }

//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

boost::shared_ptr<msTreeMapper> qcCfour::setGenBasSymLink(std::string link) {

 GenBasSymLink=link;
 return mySharedPtr();
}

//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
void qcCfour::writeInput( const msZMat& zmat) {
  
 string cmd="mkdir "+WorkingDirectory;
 struct stat st;
 if(stat(WorkingDirectory.c_str(),&st) != 0) std::system( string("mkdir "+WorkingDirectory).c_str() );

 std::system( string("(cd "+WorkingDirectory+" ; rm * )").c_str() );
 std::system( string("(cd "+WorkingDirectory+" ; ln -s "+GenBasSymLink+" )").c_str() );
/*
 if( std::system( string("mkdir "+WorkingDirectory).c_str() ) )
    BOOST_THROW_EXCEPTION( msError("can not create the working directory "+WorkingDirectory,
				  "void qcCfour::singlePoint(const msZmat& zmat)",getFullId())
			 ) ;

 if( std::system( string("(cd "+WorkingDirectory+" ; ln -s "+GenBasSymLink+" )").c_str() ) )
   BOOST_THROW_EXCEPTION( msError("can not create the symbolic link of the GENBAS file, did you set GenBasSymLink via setGenBasSymLink?",
				  "void qcCfour::singlePoint(const msZmat& zmat)",getFullId())
			 ) ;
  */ 
 int mult=1;
 if( zmat.noOfElectrons()%2 != 0) mult=2;

 string fileZmat=WorkingDirectory+"/ZMAT";
 ofstream file(fileZmat.c_str());
 
 file<<"File generated by qcCfour"<<endl;
 zmat.writeFormal(file);
 
 file<<endl<<"*CFOUR(BASIS="<<BasisSet<<",CALC="<<CiType<<",MULT="<<mult<<endl<<" CC_MAXCYC=100,LINEQ_MAXCY=100,SCF_MAXCYC=200)"<<endl<<endl;
      	//   <<"LINEQ_MAXCY=100,SCF_MAXCYC=200)"<<endl<<endl;
 
 file.close();  
}

//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

double qcCfour::singlePoint(const msZMat& zmat)
{
 writeInput(zmat);
 return singlePoint();
 }
double qcCfour::singlePoint()
{
 TestXaces2();
 string cmd ="( cd "+WorkingDirectory + "/ && xaces2 > outSinglePoint.txt)";       
 std::system(cmd.c_str());
 double e = getEnergy("outSinglePoint.txt");
 return e;
 }
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

bool  qcCfour::TestXaces2()
{
 if( std::system("xaces2") )
    BOOST_THROW_EXCEPTION( msError("can not execute xaces2; is it in your path?",
				   "bool  qcCfour::TestXaces2()",getFullId())
			  ) ;
    
 string genbas=WorkingDirectory+"/GENBAS";
 
 if( !ifstream(genbas.c_str()) )
    BOOST_THROW_EXCEPTION( msError("can not find the GENBAS file; is it in the working directory?",
				   "bool  qcCfour::TestXaces2()",getFullId())
			  ) ;
 return 1;
 }

//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

double qcCfour::getEnergy(string fileout)
{
 fileout=WorkingDirectory+"/"+fileout;
 ifstream out(fileout.c_str());
 char line[100]; string word;
  if(!out){cout<<"file "<<fileout<<" not found"<<endl; return -1;}
 while(!out.eof()) {out>>word;
                    if(word=="final")
                      {out>>word>>word>>word>>word;
	               return( atof(word.c_str()) );
	               }
                    }

 return( -1 );
 }

//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------

void qcCfour::cleanDir()
{
 
 //std::system(string(" rm "+ WorkingDirectory+"/*.*" ).c_str() );
 /* vector<string> files;
 files.push_back("/AUXMOI"); 	files.push_back("/basinfo.data"); 	files.push_back("/coriolis.joda");
 files.push_back("/DCT"); 	files.push_back("/den.dat"); 		files.push_back("/DERGAM");
 files.push_back("/DERINT");	files.push_back("/didq");		files.push_back("/DIPDER");
 files.push_back("/DIPOL"); 	files.push_back("/dnor"); 		files.push_back("/dnorraw"); 
 files.push_back("/dvr.in"); 	files.push_back("/FCM"); 		files.push_back("/FCM00"); 
 files.push_back("/FCM000"); 	files.push_back("/FCMFINAL"); 	  	files.push_back("/FCMINT"); 
 files.push_back("/FRQARC"); 	files.push_back("/FILES"); 		files.push_back("/fort.81");	
 files.push_back("/fort.82");	files.push_back("/fort.83"); 		files.push_back("/GAMLAM"); 
 files.push_back("/GRD"); 	files.push_back("/HPOLAR"); 		files.push_back("/IIII");
 files.push_back("/IIJJ");	files.push_back("/IJIJ"); 		files.push_back("/IJKL"); 
 files.push_back("/irspectrum_harmonic.txt"); 
 files.push_back("/JAINDX"); 	files.push_back("/JMOLplot"); 		files.push_back("/JOBARC"); 
 files.push_back("/MOABCD"); 	files.push_back("/MOINTS"); 		files.push_back("/MOL");  
 files.push_back("/MOLDEN");	files.push_back("/NEWMOS");   		files.push_back("/NORMCO");  
 files.push_back("/NTOTAL");  	files.push_back("/OLDMOS");   		files.push_back("/OPTARC");   	
 files.push_back("/out");  	files.push_back("/POLAR"); 		files.push_back("/QUADRATURE");  
 
 
 for(size_t i=0; i<files.size(); i++)

    if( ifstream( files[i].c_str() ) ) 
       std::system(string(" rm "+ WorkingDirectory + files[i]).c_str());  */
 }



/*
double App_C4::GetScfEnergy(string fileout)
{
  fileout=WorkDir+"/"+fileout;
 ifstream out(fileout.c_str());
 char line[100]; string word;
 if(!out){cout<<"file "<<fileout<<" not found"<<endl; return -1;}
 while(!out.eof()) {out>>word;
                    if(word=="E(SCF)=")
                      {out>>word;
	               return( atof(word.c_str()) );
	               }
                    }

 return( -1 );
 }
 */
 /*
void App_C4::OptGeom(string calcul,msZMat* zmt,string outstr,bool askFreez)
{
 if( !TestXaces2() ){ return; }
 string ans;
 vector<bool> freeze;
 freeze.resize(zmt->DoF_int(),0);
  
 if(askFreez)
  {
   std::cout<<"\033[1;34m Geometry optimization \033[0m"<<std::endl; 
   std::cout<<"Constrain some coordinates? [n]/y"<<endl;
   cin>>ans;
 
 if(ans=="y")
   {
    cout<<"enter the indice of the coordinates you want to freeze (first is 0), enter 'd' when done."<<endl;
    
    while(ans!="d"){ cin>>ans; 
	             if(ans[0]!='d'){freeze[atoi(ans.c_str())]=1;} 
	            }
    }}
    
    
 string fileZmat=WorkDir+"/ZMAT";
 std::ofstream input(fileZmat.c_str());
 input<<"File generated by Cfour package"<<endl;
 int j=-1;
 for(size_t i=0;i<zmt->Atoms.size();i++)
    {
     if(i>=0){input<<zmt->Atoms[i]->Name.c_str();}
     if(i>=1){input<<" "<<zmt->Atoms[i]->Bond.Ref+1<<" r"<<(i);j++;if( !freeze[j] ){input<<"*";}}
     if(i>=2){input<<" "<<zmt->Atoms[i]->Angle.Ref+1<<" a"<<(i-1);j++;if( !freeze[j] ){input<<"*";}}
     if(i>2){input<<" "<<zmt->Atoms[i]->Dihedre.Ref+1<<" d"<<(i-2);j++;if( !freeze[j] ){input<<"*";}}

     input<<endl;
     }
     
 input<<endl;
 for(size_t i=0;i<zmt->Atoms.size();i++)
    {if(i>=1){input<<"r"<<(i)<<" = "<<zmt->Atoms[i]->Bond.Lenght<<endl;}
     if(i>=2){input<<"a"<<(i-1)<<" = "<<zmt->Atoms[i]->Angle.Angle<<endl;}
     if(i>2){input<<"d"<<(i-2)<<" = "<<zmt->Atoms[i]->Dihedre.Angle<<endl;}
     }

 input<<endl<<"*CFOUR("<<calcul<<",MULT="<<MultStr<<endl
      <<"CmsMAXCYC=100,LINEQ_MAXCY=100,SCF_MAXCYC=200)"<<endl<<endl;
 input.close();
 
 Remove_files();
 cout<<"\033[1;34m Optimization of the structure  \033[1;0m"<<endl;
 string cmd="( cd "+WorkDir + "/ && xaces2 > "+outstr+" )";
 std::system(cmd.c_str());
 }

 
void App_C4zmat::Exe()
{ 
 ZMat=(msZMat*)Object;
 std::cout<<"Execution of the application Cfour::C4zmat"<<std::endl;
 WorkDir="molecules/"+ZMat->Parent->Name+"/data";
 string cmd="mkdir "+WorkDir;
 cmd="mkdir "+WorkDir; std::system(cmd.c_str());
 WorkDir+="/ZMat";
 cmd="mkdir "+WorkDir; std::system(cmd.c_str());
 
 string fileZmat=WorkDir+"/ZMAT"; 
 
 cout<<"Application cfour run in "<<WorkDir<<endl;
  
 if( !TestXaces2() ){ return; }
 
 //----------------------------------------
 //Create directory tree structure and file
 //----------------------------------------
 
 cout<<"Spin multiplicity?"<<endl;
 
 int mult;cin>>mult;
 MultStr="1"; if(mult==2){MultStr="2";} if(mult==3){MultStr="3";}if(mult==4){MultStr="4";}
 char a;
 cout<<"Would you like to realize a HEAT calculation ? y/[n]"<<endl;
 cin>>a;
 
 if(a=='y')
   {HEAT(); return;
    }
 cout<<"Would you like to optimize the structure ? y/[n]"<<endl;
 cin>>a;
 
 //----------------------------------------
 //Geometry optimisation
 //---------------------------------------- 
 
 if(a=='y')
   { 
    // stringstream calcul;calcul<<"*CFOUR(CALC=CCSD(T),BASIS=PVQZ,MULT="<<MultStr<<")";
     std::cout<<"Definition of the calculation type ( in parenthsesis of *CFOUR(...), without the 'mult' part)) ?"<<endl;
     string calcul;
     cin>>calcul;
     OptGeom(calcul,ZMat,"out_opt");
     UpdateZMat(ZMat);
     return;
     }
     
 //--------------------
 //Energies computation 
 //--------------------
 std::cout<<"Definition of the calculation type ( in parenthsesis of *CFOUR(...), without the 'mult' part)) ?"<<endl;
 string calcul;
 cin>>calcul;
 double E=SinglePoint(calcul,ZMat,"out");
 
 }
 
 
//------------------------------------------------------
//------------------------------------------------------
/*	    
void App_C4zmat::HEAT()
{ 
 cout<<"Would you like to optimize the structure ? y/[n]"<<endl;
 char a;cin>>a;
 
 //----------------------------------------
 //Geometry optimisation
 //---------------------------------------- 
 
 if(a=='y')
   { 
     OptGeom("CALC=CCSD(T),BASIS=PVQZ",ZMat,"out_opt");
     UpdateZMat(ZMat);

     }
     
 vector<string> calculs;vector<string> NameFile;
 calculs.push_back("CALC=HF,BASIS=AUG-PCVTZ,DBOC=ON,MULT="+MultStr); NameFile.push_back("out_HF3");
 calculs.push_back("CALC=CCSD(T),BASIS=AUG-PCVTZ,MULT="+MultStr); NameFile.push_back("out_CCSDt_3");
 calculs.push_back("CALC=CCSD(T),BASIS=AUG-PCVTZ,RELATIVITSIC=MVD1,MULT="+MultStr); NameFile.push_back("out_CCSDt_3_REL");
 calculs.push_back("CALC=CCSD(T),BASIS=AUG-PCVQZ,MULT="+MultStr); NameFile.push_back("out_CCSDt_4");
 calculs.push_back("CALC=CCSD(T),BASIS=AUG-PCV5Z,MULT="+MultStr);NameFile.push_back("out_CCSDt_5"); 
 calculs.push_back("CALC=CCSD(T),BASIS=PVTZ,FROZEN_CORE=ON,MULT="+MultStr);NameFile.push_back("out_CCSDt_3_FC");
 calculs.push_back("CALC=CCSD(T),BASIS=PVQZ,FROZEN_CORE=ON,MULT="+MultStr);NameFile.push_back("out_CCSDt_4_FC");
 calculs.push_back("CALC=CCSDT,BASIS=PVTZ,FROZEN_CORE=ON,MULT="+MultStr);NameFile.push_back("out_CCSDT_3_FC");
 calculs.push_back("CALC=CCSDT,BASIS=PVQZ,FROZEN_CORE=ON,MULT="+MultStr);NameFile.push_back("out_CCSDT_4_FC");
 calculs.push_back("CALC=CCSDT,BASIS=PVDZ,FROZEN_CORE=ON,MULT="+MultStr);NameFile.push_back("out_CCSDT_2_FC");  
 calculs.push_back("CALC=CCSDT(Q),BASIS=PVDZ,FROZEN_CORE=ON,MULT="+MultStr);NameFile.push_back("out_CCSDTq_2_FC");  
  
 for(size_t i=0;i<calculs.size();i++)
    {
     string filename=WorkDir+"/"+NameFile[i];
     if(  ifstream(filename.c_str()) && (GetEnergy(NameFile[i])!=-1) )
       {
        cout<<"\033[1;34m Calculation "<<i<<"/"<<calculs.size()<<": "<<calculs[i]<<" already done, found the file with the energy: "<<GetEnergy(NameFile[i])<<"\033[1;0m"<<endl;
        
        }
     else{
          string fileZmat=WorkDir+"/ZMAT";
          ofstream file(fileZmat.c_str());
          file<<"File generated by App_C4Scan"<<endl;
          ZMat->WriteFormal(file);
          file<<endl<<"*CFOUR("<<calculs[i]<<")"<<endl<<endl;
          file.close();
          Remove_files();
          cout<<"\033[1;34m Calculation "<<i+1<<"/"<<calculs.size()<<":"<<calculs[i]<<" \033[1;0m"<<endl;
          string cmd="( cd "+WorkDir + "/ && xaces2 > "+NameFile[i]+" )";
          std::system(cmd.c_str());
	  double e=GetEnergy(NameFile[i]);
	  cout<<"\033[1;32m \t electronic energy= "<<e<<" hartree/part \033[1;0m"<<endl;
          }
     }
 Remove_files();
 cout<<endl<<setprecision(8)<<"calculation done"<<endl;
 


 msheatEnergy HEAT;
 
 msccFellerExtrapol* E_HF=new msccFellerExtrapol();
 E_HF->i0=3;
 E_HF->e0=GetEnergy("out_HF3");
 E_HF->e1=GetScfEnergy("out_CCSDt_4");
 E_HF->e2=GetScfEnergy("out_CCSDt_5"); 
 HEAT.E_HF=E_HF;
 
 cout<<"E_HF inf="<<HEAT.E_HF->getExtrapolValue()<<endl;
 
 msccApweExtrapol* Delta_CCSDt=new msccApweExtrapol();
 Delta_CCSDt->i0=4;
 Delta_CCSDt->e0=GetEnergy("out_CCSDt_4")-GetScfEnergy("out_CCSDt_4");
 Delta_CCSDt->e1=GetEnergy("out_CCSDt_5")-GetScfEnergy("out_CCSDt_5"); 
 HEAT.Delta_CCSDt=Delta_CCSDt;
 
 cout<<"DeltaE_CCSDt inf="<<HEAT.Delta_CCSDt->getExtrapolValue()<<endl;
  
 msccApweExtrapol* E_CCSDT=new msccApweExtrapol();
 E_CCSDT->i0=3;
 E_CCSDT->e0=GetEnergy("out_CCSDT_3_FC");
 E_CCSDT->e1=GetEnergy("out_CCSDT_4_FC"); 
 HEAT.E_CCSDT=E_CCSDT;
 
 msccApweExtrapol* E_CCSDt=new msccApweExtrapol();
 E_CCSDt->i0=3;
 E_CCSDt->e0=GetEnergy("out_CCSDt_3_FC");
 E_CCSDt->e1=GetEnergy("out_CCSDt_4_FC"); 
 HEAT.E_CCSDt=E_CCSDt;
 
 cout<<"E_CCSDT-E_CCSDt inf="<<HEAT.E_CCSDT->getExtrapolValue()-HEAT.E_CCSDt->getExtrapolValue()<<endl;

 cout<<"E_CCSDTQ-E_CCSDT inf="<<GetEnergy("out_CCSDTq_2_FC")-GetEnergy("out_CCSDT_2_FC")<<endl;
 


 cout<<"Xml node to include in the species definition, DBOC and SO has to be included by hand form the files out_CCSDt_3_REL and out_HF3:"<<endl<<endl;
 cout<<setprecision(8)<<endl;
 cout<<" <heatEnergy-E0s DBOC=\"0\" REL=\"0\" > "<<endl;
 cout<<"    <Feller-HF i0=\"3\" e0=\""<< E_HF->e0<<"\" e1=\""<< E_HF->e1<<"\" e2=\""<<E_HF->e2<<"\" Comment=\"BASIS=AUG-PCVxZ\" /> "<<endl;
 cout<<"    <Apwe-delta_CCSDt i0=\"4\" e0=\""<< Delta_CCSDt->e0<<"\" e1=\""<< Delta_CCSDt->e1<<"\"  Comment=\"BASIS=AUG-PCVxZ\" /> "<<endl;
 cout<<"    <Apwe-CCSDT i0=\"3\" e0=\""<< E_CCSDT->e0<<"\" e1=\""<< E_CCSDT->e1<<"\" Comment=\"BASIS=PVxZ and frozen core\"/> "<<endl;
 cout<<"    <Apwe-CCSDt i0=\"3\" e0=\""<<  E_CCSDt->e0<<"\" e1=\""<< E_CCSDt->e1<<"\" Comment=\"BASIS=PVxZ and frozen core\" /> "<<endl;
 cout<<"    <Single-HLC i0=\"2\" e0=\""<< GetEnergy("out_CCSDTq_2_FC")-GetEnergy("out_CCSDT_2_FC") <<"\" Comment=\"CCSDT(Q)//PVDZ and frozen core\"/> "<<endl;
 cout<<" </heatEnergy-E0s> "<<endl;  
 }
*/



}
